#tablice wielowymiarowe w Julii
Asmall=[[1.0 0.0 10]; [0.0 1.0 10]]
Bsmall=Asmall
Asmall
C = zeros(Float64, size(Asmall,1), size(Asmall,2))

size(Asmall)

# mnożenie macierzy - wersja naiwna, naiwna przez sprosob dostepu do pamieci
# najlepiej zeby dane byly blisko siebie w pamieci, przez sposob dzialania pamieci cache, ktora odczytuje dane 
# blokami zwanymi liniami cache. Odczytywanie obok siebie jest duzo szybsze.
# Pytniem jest jak jest tablica przechowywana w pamieci, zazwyczaj wierszami (C) ale w julli sa przechowywane
# kolumnami 
function naive_multiplication(A,B)
C=zeros(Float64,size(A,1),size(B,2))
  for i=1:size(A,1)
    for j=1:size(B,2)
        for k=1:size(A,2)
            C[i,j]=C[i,j]+A[i,k]*B[k,j]
        end
    end
end
C
end

#kompilacja
naive_multiplication(Asmall,Bsmall)

#kompilacja funkcji BLASowej do mnożenia macierzy
#https://docs.julialang.org/en/stable/stdlib/linalg/#BLAS-Functions-1
#to tez jest mnozenie macierzy ale zoptymalizowanymi funkcjami BLAS
Asmall*Bsmall

A=rand(1000,1000); #tworzenie macierzy 1000 x 1000 z losowymi wartosciami 
B=rand(1000,1000);

# Należy pamiętać o "column-major" dostępie do tablic - 
# pierwszy indeks zmienia się szybciej
# tak jak Matlab, R, Fortran 
# inaczej niz C, Python
A1 = [[1 2]; [3 4]]
vec(A1)

# poprawiona funkcja korzytająca z powyższego oraz z faktu, że
#można zmieniać kolejność operacji dodawania (a co za tym idzie kolejnosc petli).
# jest lepsza przez zamienienie kolejnosci petli i czesciej odczytuje elementy bedace blizej siebie 
function better_multiplication( A,B )
C=zeros(Float64,size(A,1),size(B,2))
  for j=1:size(B,2)
    for k=1:size(A,2)
        for i=1:size(A,1)
            C[i,j]=C[i,j]+A[i,k]*B[k,j]
        end
    end
end
C
end

better_multiplication(Asmall, Bsmall)

@elapsed naive_multiplication(A,B) #mierzenie czasu

@elapsed better_multiplication(A,B)

@elapsed A*B #blas

# aproksymacja sredniokwadratowa wielomianem - tutaj przyklad dla wielomianu 3 stopnia
# pakiet Polynomials jest mozliwy do instalacji pod Juliabox
# https://github.com/JuliaMath/Polynomials.jl


using Polynomials
xs = 0:3; ys = [1,3,4,5]
fit1=polyfit(xs, ys,3)

# po prostu za xs podstawic to co mam za wielkosc macierzy a za ys podstawiac po kolei te zmienne dla których
# chce wyliczyc wielomian

# obliczanie wartosci wielomianu 
fit1(1)

# obliczanie wartosci wielomianu (drugi sposób)
polyval(fit1, 1)

using Plots

# geste punkty do wyliczenia wartosci wielomianu aproksymujacego:
xd=0:0.1:10
# wykres wartosci wielomianu dla gestych punktow:
plot(xd,polyval(fit1, xd))

# ! -dodanie do tego samego wykresu punktów wg ktorych aproksymowalismy
scatter!(xs,ys)

columns_and_rows = Int64[]
naive_time = Float64[]
better_time = Float64[]
blas_time = Float64[]

nb_of_tests = 10
i = 50
while(i <= 1000)
    for k=0:nb_of_tests
        A = rand(i,i)
        B = rand(i,i)
        push!(columns_and_rows, i)
        push!(naive_time,@elapsed naive_multiplication(A,B))
        push!(better_time,@elapsed better_multiplication(A,B))
        push!(blas_time,@elapsed A * B)
    end
    i += 50
end
columns_and_rows

using DataFrames

df = DataFrame()
df[:columns_and_rows]= columns_and_rows
df[:naive_time] = naive_time
df[:better_time] = better_time
df[:blas_time] = blas_time

using Statistics
df2 = DataFrame(by(df, [:columns_and_rows],
    :naive_time => mean,
    :naive_time => std,
    :better_time => mean,
    :better_time => std,
    :blas_time => mean,   
    :blas_time => std))

using Plots

ydata = scatter(df2[:columns_and_rows],
    [df2[:naive_time_mean] df2[:better_time_mean] df2[:blas_time_mean]], 
    yerr = [df2[:naive_time_std] df2[:better_time_std] df2[:blas_time_std]],
    labels = ["Naive Time" "Better time" "BLAS time"],
    title = "Time of mulitplication in Julia",
    legend=:topleft,
    xlabel = "Size of matrix",
    ylabel = "Time [s]",)

using CSV
input0="resultO0.csv"
mydata0=CSV.read(input0, delim=";")
input1="resultO1.csv"
mydata1=CSV.read(input1, delim=";")
input2="resultO2.csv"
mydata2=CSV.read(input2, delim=";")
input3="resultO3.csv"
mydata3=CSV.read(input3, delim=";")

using Statistics, DataFrames, Plots
data0 = DataFrame(by(mydata0, [:columns_and_rows],
    :naive_time => mean,
    :naive_time => std,
    :better_time => mean,
    :better_time => std,
    :blas_time => mean,   
    :blas_time => std))
data1 = DataFrame(by(mydata1, [:columns_and_rows],
    :naive_time => mean,
    :naive_time => std,
    :better_time => mean,
    :better_time => std,
    :blas_time => mean,   
    :blas_time => std))
data2 = DataFrame(by(mydata2, [:columns_and_rows],
    :naive_time => mean,
    :naive_time => std,
    :better_time => mean,
    :better_time => std,
    :blas_time => mean,   
    :blas_time => std))
data3 = DataFrame(by(mydata3, [:columns_and_rows],
    :naive_time => mean,
    :naive_time => std,
    :better_time => mean,
    :better_time => std,
    :blas_time => mean,   
    :blas_time => std))

plotO0 = scatter(data0[:columns_and_rows],
    [data0[:naive_time_mean] data0[:better_time_mean] data0[:blas_time_mean]], 
    yerr = [data0[:naive_time_std] data0[:better_time_std] data0[:blas_time_std]],
    labels = ["Naive Time" "Better time" "BLAS time"],
    title = "Time of mulitplication in C, with O0 optim",
    legend=:topleft,
    xlabel = "Size of matrix",
    ylabel = "Time [s]")
plotO1 = scatter(data1[:columns_and_rows],
    [data1[:naive_time_mean] data1[:better_time_mean] data1[:blas_time_mean]], 
    yerr = [data1[:naive_time_std] data1[:better_time_std] data1[:blas_time_std]],
    labels = ["Naive Time" "Better time" "BLAS time"],
    title = "Time of mulitplication in C, with O1 optim",
    legend=:topleft,
    xlabel = "Size of matrix",
    ylabel = "Time [s]")
plotO2 = scatter(data2[:columns_and_rows],
    [data2[:naive_time_mean] data2[:better_time_mean] data2[:blas_time_mean]], 
    yerr = [data2[:naive_time_std] data2[:better_time_std] data2[:blas_time_std]],
    labels = ["Naive Time" "Better time" "BLAS time"],
    title = "Time of mulitplication in C, with O2 optim",
    legend=:topleft,
    xlabel = "Size of matrix",
    ylabel = "Time [s]")

plotO3 = scatter(data3[:columns_and_rows],
    [data3[:naive_time_mean] data3[:better_time_mean] data3[:blas_time_mean]], 
    yerr = [data3[:naive_time_std] data0[:better_time_std] data3[:blas_time_std]],
    labels = ["Naive Time" "Better time" "BLAS time"],
    title = "Time of mulitplication in C, with O3 optim",
    legend=:topleft,
    xlabel = "Size of matrix",
    ylabel = "Time [s]")


naive_compare = scatter([data0[:columns_and_rows] data1[:columns_and_rows] data2[:columns_and_rows]],
    [data0[:naive_time_mean] 
     data1[:naive_time_mean] 
     data2[:naive_time_mean] ],
    yerr = 
    [data0[:naive_time_std] 
     data1[:naive_time_std] 
     data2[:naive_time_std] ],
    labels = ["Naive Time O0" "Naive Time O1" "Naive Time O2"],
    title = "naive in optimalization dependency",
    legend=:topleft,
    xlabel = "Size of matrix",
    ylabel = "Time [s]",)


display(plotO0)
display(plotO1)
display(plotO2)
display(plotO3)

display(naive_compare)

using Polynomials, Plots

fit_data0_naive=polyfit(data0[:columns_and_rows],data0[:naive_time_mean],3)
xd=0:0.01:1000
plot(xd,polyval(fit_data0_naive, xd))
scatter!(data0[:columns_and_rows],data0[:naive_time_mean])
scatter!(labels = ["Naive Time O0" "Polynomial approximation"],
legend=:topleft,xlabel = "Size of matrix", ylabel = "Time [s]", title = "Naive in C")
#print(fit_data0_naive)

fit_data0_better=polyfit(data0[:columns_and_rows],data0[:better_time_mean],3)
xd=0:0.01:1000
plot(xd,polyval(fit_data0_better, xd))
scatter!(data0[:columns_and_rows],data0[:naive_time_mean])
scatter!(labels = ["Naive Time O0" "Polynomial approximation"],
legend=:topleft,xlabel = "Size of matrix", ylabel = "Time [s]", title = "Better in C")
#print(fit_data0_better)

fit_data0_blas=polyfit(data0[:columns_and_rows],data0[:blas_time_mean],3)
xd=0:0.01:1000
plot(xd,polyval(fit_data0_blas, xd))
scatter!(data0[:columns_and_rows],data0[:naive_time_mean])
scatter!(labels = ["Naive Time O0" "Polynomial approximation"],
legend=:topleft,xlabel = "Size of matrix", ylabel = "Time [s]", title = "BLAS in C")
#print(fit_data0_blas)

fit_julia_naive=polyfit(df2[:columns_and_rows],df2[:naive_time_mean],3)
xd=0:0.01:1000
plot(xd,polyval(fit_data0_naive, xd))
scatter!(data0[:columns_and_rows],data0[:naive_time_mean])
scatter!(labels = ["Naive Time O0" "Polynomial approximation"],
legend=:topleft,xlabel = "Size of matrix", ylabel = "Time [s]", title = "Naive in Julia")
#print(fit_julia_naive)

fit_julia_better=polyfit(df2[:columns_and_rows],df2[:better_time_mean],3)
xd=0:0.01:1000
plot(xd,polyval(fit_data0_better, xd))
scatter!(data0[:columns_and_rows],data0[:naive_time_mean])
scatter!(labels = ["Naive Time O0" "Polynomial approximation"],
legend=:topleft,xlabel = "Size of matrix", ylabel = "Time [s]", title = "Better in Julia")
#print(fit_julia_better)

fit_julia_blas=polyfit(df2[:columns_and_rows],df2[:blas_time_mean],3)
xd=0:0.01:1000
plot(xd,polyval(fit_data0_blas, xd))
scatter!(data0[:columns_and_rows],data0[:naive_time_mean])
scatter!(labels = ["Naive Time O0" "Polynomial approximation"],
legend=:topleft,xlabel = "Size of matrix", ylabel = "Time [s]", title = "BLAS in Julia")
#print(fit_julia_blas)